#include <navierstokes.hpp>
#include <life/lifefilters/geotool.hpp>
#include <iostream>

namespace Life {

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
typename NavierStokes<Dim,Order,Entity>::mesh_ptrtype
NavierStokes<Dim,Order,Entity>::createMesh( double meshSize )
{

    mesh_ptrtype mesh( new mesh_type );

    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->application()->vm()["hsize"].template as<double>()
                                           );


    //Geometry
    GeoTool::Node x1(0,0);
    GeoTool::Node x2(4,1);
    GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
    R.setMarker(_type="line",_name="Entree",_marker4=true);
    //R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker3=true);
    R.setMarker(_type="line",_name="ParoiB",_marker1=true);
    R.setMarker(_type="line",_name="ParoiH",_marker3=true);
    R.setMarker(_type="line",_name="Sortie",_marker2=true);
    R.setMarker(_type="surface",_name="OmegaFluide",_marker1=true);

    GeoTool::Node x_centre(1,0.5);
    GeoTool::Node x_bord(1.1,0.5);
    GeoTool::Circle C( meshSize/6.,"disque",x_centre,x_bord);
    C.setMarker(_type="line",_name="Paroi",_markerAll=true);
    C.setMarker(_type="surface",_name="OmegaDisque",_marker1=true);

    //std::cout<< (R-C).geoStr();

    Gmsh gmsh;
    std::string fname = gmsh.generate( "domain", (R-C).geoStr()  );

    ImporterGmsh<mesh_type> import( fname );
    import.setVersion( "2.1" );
    mesh->accept( import );

    return mesh;

} // Stokes::createMesh


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
typename NavierStokes<Dim,Order,Entity>::mesh_structure_ptrtype
NavierStokes<Dim,Order,Entity>::createStructMesh( double meshSize )
{

  mesh_structure_ptrtype struct_mesh( new mesh_structure_type );
  GmshTensorizedDomain<1,1,1,Simplex> td;
  td.setCharacteristicLength( meshSize );
  td.setX( std::make_pair( 0, 4 ) );
  std::string fname = td.generate( "boundary_description" );
  ImporterGmsh<mesh_structure_type> struct_import( fname );
  struct_mesh->accept( struct_import );
  return struct_mesh;

}



template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::solve( sparse_matrix_ptrtype const& D,
                                   element_fluid_type& u,
                                   vector_ptrtype const& F,
                                   bool is_sym )
{
    vector_ptrtype U( M_backend->newVector( u.functionSpace() ) );
    M_backend->solve( D, D, U, F, false );
    u = *U;
} // NavierStokes::solve

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::exportResults( element_fluid_type& U )
{
    exporter->step( 0 )->setMesh( U.functionSpace()->mesh() );
    exporter->step( 0 )->add( "u", U.template element<0>() );
    exporter->step( 0 )->add( "p", U.template element<1>() );
    exporter->save();
} // NavierStokes::export


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::exportResultsInstat( element_fluid_type& U,element_fluid_velocity_type & mesh_velocity,
						       ale_map_element_type & aux,double time )
{

  this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/Export/" )
			  % this->application()->about().appName()
			  % convex_type::name()
			  % Order
			  % this->application()->vm()["hsize"].template as<double>()
			  );


  //exporter->step( time )->setMesh( U.functionSpace()->mesh() );
  //exporter->step( time )->setMesh( mesh_velocity.functionSpace()->mesh() );
  exporter->step( time )->setMesh( mesh_fluid );
  exporter->step( time )->add( "u", U.template element<0>() );
  exporter->step( time )->add( "p", U.template element<1>() );
  exporter->step( time )->add( "mesh_velocity", mesh_velocity );
  exporter->step( time )->add( "aux", aux );
  exporter->save();
  ++M_cpt_export;
} // NavierStokes::export

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
NavierStokes<Dim, Order, Entity>::exportResultsInstat( mesh_ptrtype __mesh,element_fluid_type& U,element_fluid_velocity_type & mesh_velocity,
						       ale_map_element_type & aux,double time )
{

  this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/Export/" )
                                         % this->application()->about().appName()
                                         % convex_type::name()
                                         % Order
                                         % this->application()->vm()["hsize"].template as<double>()
                                         );

  exporter->step( time )->setMesh( __mesh );
  exporter->step( time )->add( "u", U.template element<0>() );
  exporter->step( time )->add( "p", U.template element<1>() );
  exporter->step( time )->add( "mesh_velocity", mesh_velocity );
  exporter->step( time )->add( "aux", aux );
  exporter->save();

  ++M_cpt_export;
} // NavierStokes::export

template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
double
NavierStokes<Dim, Order, Entity>::diffL2( element_fluid_type& U1, element_fluid_type& U2 )
{
  using namespace Life::vf;

  element_fluid_0_type u1 = U1.template element<0>();
  element_fluid_1_type p1 = U1.template element<1>();
  element_fluid_0_type u2 = U2.template element<0>();
  element_fluid_1_type p2 = U2.template element<1>();

  mesh_ptrtype mesh = U1.functionSpace()->mesh();

  double __norm = math::sqrt( integrate( elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
					 trans(idv(u1))*idv(u1) +idv(p1)*idv(p1) ).evaluate()(0,0) );
  double __diff = math::sqrt( integrate(elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
					trans(idv(u1)-idv(u2))*(idv(u1)-idv(u2))
					+(idv(p1)-idv(p2))*(idv(p1)-idv(p2)) ).evaluate()(0,0) );
  if ( __norm != 0 ) __diff /= __norm;

  return __diff;
}

} // Life

// instantiation
template class Life::NavierStokes<2,2,Life::Simplex>;
//template class Life::NavierStokes<3,2,Life::Simplex>;

